Specifically, we will explore the results of 827 datapoints from 28 studies conducted in 3 oceans on 87 species within 46 genera.
## n
## 1 827
## Ref n
## 1 DS08a 162
## 2 DS08b 113
## 3 DS71 112
## 4 DS04 70
## 5 DS68 48
## 6 DS46 44
## 7 DS29 40
## 8 DS30 36
## 9 DS17 24
## 10 DS45 24
## 11 DS16 19
## 12 DS18 16
## 13 DS27 16
## 14 DS32 16
## 15 DS07 13
## 16 DS03 12
## 17 DS43 12
## 18 DS05 10
## 19 DS33 8
## 20 DS12 7
## 21 DS26 6
## 22 DS11 4
## 23 DS21 4
## 24 DS28 4
## 25 DS25 2
## 26 DS36 2
## 27 DS69 2
## 28 DS23 1
## Ref_name n
## 1 Hodgson (1989) PhD Chapter IV 275
## 2 HDR EOC and CSA Ocean Services, Inc. (2014) 112
## 3 Duckworth et al. (2017) 70
## 4 Flores et al. (2012) 48
## 5 Stafford-Smith (1993) 44
## 6 Stafford-Smith (1990) PhD Chapter 4 40
## 7 Stafford-Smith (1992) 36
## 8 Piniak (2007) 24
## 9 Rogers (1983) 24
## 10 Philipp and Fabricius (2003) 19
## 11 Piniak and Brown (2008) 16
## 12 Sofonia (2006) PhD Chapter 4 16
## 13 Stewart et al. (2006) 16
## 14 Hodgson (1990b) 13
## 15 Bessell-Browne et al. (2017a) 12
## 16 Peters and Pilson (1985) 12
## 17 Hodel (2007) 10
## 18 Vargas-Angel et al. (2006) 8
## 19 Loiola et al. (2013) 7
## 20 Sofonia (2006) PhD Chapter 3 6
## 21 Lirman et al. (2008) 4
## 22 Rogers (1979) 4
## 23 Sofonia and Anthony (2008) 4
## 24 Gil et al. (2016) 2
## 25 Shore-Maggio et al. (2018) 2
## 26 Zill et al. (2017) 2
## 27 Selim (2007) 1
## Ocean n
## 1 Pacific 757
## 2 Atlantic 69
## 3 Indian 1
## Selecting by n
## Gsp n
## 1 Acropora_millepora 38
## 2 Leptoria_phrygia 42
## 3 Montipora_aequituberculata 26
## 4 Oxypora_glabra 23
## 5 Pavona_cactus 29
## 6 Porites_cylindrica 30
## 7 Porites_lobata/lutea 56
## 8 Porites_lutea 28
## 9 Porites_rus 30
## 10 Turbinaria_reniformis 24
## Selecting by n
## Updated_Genus n
## 1 Porites 152
## 2 Montipora 137
## 3 Acropora 88
## 4 Leptoria 42
## 5 Pavona 37
## 6 Turbinaria 33
## 7 Oxypora 31
## 8 Pocillopora 29
## 9 Echinopora 23
## 10 Mycedium 21
Specifically, we will explore the results of 272 datapoints from 6 studies conducted in 3 oceans on 19 species within 14 genera (shown below only the top 10 species and genera, by number of datapoints).
## n
## 1 272
## Ref n
## 1 SS06 150
## 2 SS08 48
## 3 SS22c 28
## 4 SS27 16
## 5 SS22a 10
## 6 SS04 9
## 7 SS22b 8
## 8 SS25 3
## Ref_name n
## 1 Browne et al. 2015 150
## 2 Flores et al. 2012 48
## 3 Rice 1984 46
## 4 Anthony et al. 2007 16
## 5 Bessell-Browne et al. 2017 9
## 6 Te 2001 Chapter 6 3
## Ocean n
## 1 Indian 150
## 2 Pacific 76
## 3 Atlantic 46
## Selecting by n
## Gsp n
## 1 Merulina_ampliata 50
## 2 Pachyseris_speciosa 50
## 3 Platygyra_sinensis 50
## 4 Montipora_aequituberculata 30
## 5 Acropora_millepora 21
## 6 Acropora_intermedia 16
## 7 Cladocora_arbuscula 8
## 8 Manicina_areolata 8
## 9 Solenastrea_hyades 8
## 10 Siderastrea_radians 6
## Selecting by n
## Updated_Genus n
## 1 Merulina 50
## 2 Pachyseris 50
## 3 Platygyra 50
## 4 Acropora 37
## 5 Montipora 36
## 6 Cladocora 8
## 7 Manicina 8
## 8 Solenastrea 8
## 9 Siderastrea 6
## 10 Isophyllia 4
## 11 Phyllangia 4
## 12 Scolymia 4
## 13 Stephanocoenia 4
First let’s explore all data from all species for which there is binary data about the presence of ‘small necroses’, ‘large necroses’, and ‘total mortality’ as a result of exposure to deposited sediment.
Now let’s calculate the thresholds for each category of mortality, based on the binary data explored above.
## LOAEL
## 1 4.9
## NOAEL
## 1 4.4
## LOAEL
## 1 0.917
## NOAEL
## 1 0.917
## LOAEL
## 1 20.8
## NOAEL
## 1 20.8
## LOAEL
## 1 3
## NOAEL
## 1 3
## LOAEL
## 1 20.8
## NOAEL
## 1 20.8
## LOAEL
## 1 1
## NOAEL
## 1 1
## LOAEL
## 1 3.22
## NOAEL
## 1 3.22
## LOAEL
## 1 14
## NOAEL
## 1 14
## LOAEL
## 1 29.1
## NOAEL
## 1 29.1
## LOAEL
## 1 84
## NOAEL
## 1 84
## LOAEL
## 1 29.1
## NOAEL
## 1 29.1
## LOAEL
## 1 84
## NOAEL
## 1 84
## n
## 1 602
## Data: smallDS2
## Models:
## modDS6: Binary_small_necroses ~ Sed_level_stand_mg * Sed_exposure_days +
## modDS6: (1 | Gsp)
## modDS9: Binary_small_necroses ~ Sed_level_stand_mg * Sed_exposure_days +
## modDS9: (1 | Ref)
## modDS12: Binary_small_necroses ~ Sed_level_stand_mg * Sed_exposure_days +
## modDS12: (1 | Gsp) + (1 | Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modDS6 5 660.76 682.77 -325.38 650.76
## modDS9 5 702.81 724.81 -346.40 692.81 0.000 0 1
## modDS12 6 652.31 678.71 -320.15 640.31 52.496 1 4.31e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Data: smallDS2
## Models:
## modDS10: Binary_small_necroses ~ Sed_level_stand_mg + (1 | Gsp) + (1 |
## modDS10: Ref)
## modDS11: Binary_small_necroses ~ Sed_level_stand_mg + Sed_exposure_days +
## modDS11: (1 | Gsp) + (1 | Ref)
## modDS12: Binary_small_necroses ~ Sed_level_stand_mg * Sed_exposure_days +
## modDS12: (1 | Gsp) + (1 | Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modDS10 4 657.21 674.81 -324.61 649.21
## modDS11 5 650.97 672.97 -320.48 640.97 8.2432 1 0.00409 **
## modDS12 6 652.31 678.71 -320.15 640.31 0.6604 1 0.41641
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Binary_small_necroses ~ Sed_level_stand_mg + Sed_exposure_days +
## (1 | Gsp) + (1 | Ref)
## Data: smallDS2
##
## AIC BIC logLik deviance df.resid
## 651.0 673.0 -320.5 641.0 597
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.1315 -0.6021 -0.1801 0.5730 3.0137
##
## Random effects:
## Groups Name Variance Std.Dev.
## Gsp (Intercept) 2.871 1.694
## Ref (Intercept) 2.404 1.551
## Number of obs: 602, groups: Gsp, 75; Ref, 21
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.975421 0.562799 -3.510 0.000448 ***
## Sed_level_stand_mg 0.002889 0.001178 2.453 0.014167 *
## Sed_exposure_days 0.013832 0.004908 2.818 0.004828 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Binary_small_necroses ~ log10_conc + log10_dur + (1 | Gsp) +
## (1 | Ref)
## Data: smallDS2
##
## AIC BIC logLik deviance df.resid
## 639.8 661.8 -314.9 629.8 597
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.3280 -0.5459 -0.1809 0.5600 3.1754
##
## Random effects:
## Groups Name Variance Std.Dev.
## Gsp (Intercept) 3.219 1.794
## Ref (Intercept) 3.934 1.983
## Number of obs: 602, groups: Gsp, 75; Ref, 21
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.4867 1.1307 -4.852 1.22e-06 ***
## log10_conc 1.3053 0.3869 3.374 0.000742 ***
## log10_dur 1.5651 0.4175 3.749 0.000178 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] 0.8272425
##
## p 0 1
## 0 344 61
## 1 43 154
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
## Area under the curve: 0.9049
## R2m R2c
## theoretical 0.07964959 0.7100653
## delta 0.07245285 0.6459074
## $Gsp
##
## $Ref
## Groups Name Std.Dev.
## Gsp (Intercept) 1.7942
## Ref (Intercept) 1.9834
## [1] 7.267637
## Est LL UL
## (Intercept) 0.004141512 0.0004515197 0.03798754
## log10_conc 3.688788053 1.7279975250 7.87452360
## log10_dur 4.783110717 2.1103679520 10.84083376
For every 10-fold increase in exposure concentration of deposited sediment, the odds of developing small necroses (<50% of coral tissue surface area) increase by 3.7 times (95% CI 1.7, 7.9; GLMM z = 3.374, p = 0.0007), after accounting for exposure duration and variability among studies and species. For every 10-fold increase in exposure duration of deposited sediment, the odds of developing small necroses increase by 4.8 times (95% CI 2.1, 10.8; GLMM z = 3.749, p = 0.0002), after accounting for exposure concentration and variability among studies and species.
#Overlaying glmm results on figure, but prediction line is jagged...
pred_modDS11_log <- predict(modDS11_log, newdata=smallDS2, type="response")
smallDS3 <- cbind(smallDS2, pred_modDS11_log)
smallDS_plot <- ggplot(data = smallDS3) +
geom_point(mapping = aes(
x = Sed_level_stand_mg,
y = Binary_small_necroses)) +
labs(x = expression("Sediment exposure concentration (mg/cm"^"2"*"/day)"),
y = "Small necroses due to sediment exposure?") +
scale_x_log10(limits=c(1,1000),breaks=c(0.01,0.1,1,10,100,1000),
label=c("0.01","0.1","1","10","100","1000")) +
geom_line(aes(x = Sed_level_stand_mg, y = pred_modDS11_log), inherit.aes=FALSE) +
theme_bw()
smallDS_plot
That plot is confusing to interpret. Let’s see if I can plot the average marginal probability, i.e., the average change in probability of the outcome across the range of the predictor of interest. This is described in some detail at the following, useful website: https://stats.idre.ucla.edu/r/dae/mixed-effects-logistic-regression/
jvalues <- with(smallDS2, seq(from = min(log10_conc), to = max(log10_conc), length.out = 100))
# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
smallDS2$log10_conc <- j
predict(modDS11_log, newdata = smallDS2, type = "response")
})
# average marginal predicted probability across a few different sediment levels
#sapply(pp[c(1, 20, 40, 60, 80, 100)], mean)
##
# get the means with lower and upper quartiles
plotdat <- t(sapply(pp, function(x) {
c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))
# add in log10_conc values and convert to data frame
plotdat <- as.data.frame(cbind(plotdat, jvalues))
plotdat <- plotdat %>% mutate(Sed_level_linear=10^jvalues)
# better names and show the first few rows
colnames(plotdat) <- c("PredictedProbability", "MedianProbability",
"LowerQuantile", "UpperQuantile", "log10_conc", "Sed_conc")
#head(plotdat)
# plot average marginal predicted probabilities
ggplot(plotdat, aes(x = log10_conc)) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
ggplot(plotdat, aes(x = log10_conc)) +
geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
ggplot(plotdat, aes(x = Sed_conc, y = PredictedProbability)) +
geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1)) +
scale_x_log10()
#Overlaying average marginal predicted probabilities on figure
#by Ref
smallDS_plot2 <- ggplot() +
geom_point(data = smallDS2, mapping = aes(
x = Sed_level_stand_mg,
y = Binary_small_necroses,
color = Ref)) +
labs(x = expression("Sediment exposure concentration (mg/cm"^"2"*"/day)"),
y = "Predicted probability of small\nnecroses due to sediment exposure",
color = "\nBinary Data by Study",
linetype = "Predicted\nProbability") +
scale_x_log10(limits=c(1,1000), breaks=c(0.01,0.1,1,10,100,1000,loaelDS1),
label=c("0.01","0.1","1","10","100","1000",round(loaelDS1,digits=1))) +
geom_ribbon(data = plotdat, aes(x = Sed_conc, y = PredictedProbability,
ymin = LowerQuantile, ymax = UpperQuantile),
alpha = 0.15) +
geom_line(data = plotdat, aes(x = Sed_conc, y = PredictedProbability,
linetype = "twodash")) +
geom_line(data = plotdat, aes(x = Sed_conc, y = MedianProbability,
linetype = "solid")) +
geom_vline(xintercept=loaelDS1, linetype="dashed", color = "red") +
theme_bw() +
scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))
smallDS_plot2
jvalues <- with(smallDS2, seq(from = min(log10_dur), to = max(log10_dur), length.out = 100))
# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
smallDS2$log10_dur <- j
predict(modDS11_log, newdata = smallDS2, type = "response")
})
# average marginal predicted probability across a few different sediment levels
#sapply(pp[c(1, 20, 40, 60, 80, 100)], mean)
## [1] 0.1797186 0.1960244 0.2142042 0.2333993 0.2535694 0.2746604
# get the means with lower and upper quartiles
plotdat <- t(sapply(pp, function(x) {
c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))
# add in Sed_exposure_days values and convert to data frame
plotdat <- as.data.frame(cbind(plotdat, jvalues))
plotdat <- plotdat %>% mutate(Sed_dur_linear=10^jvalues)
# better names and show the first few rows
colnames(plotdat) <- c("PredictedProbability", "MedianProbability",
"LowerQuantile", "UpperQuantile", "log10_dur","Sed_dur")
#head(plotdat)
# plot average marginal predicted probabilities
ggplot(plotdat, aes(x = log10_dur)) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
ggplot(plotdat, aes(x = log10_dur)) +
geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
ggplot(plotdat, aes(x = Sed_dur, y = PredictedProbability)) +
geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1)) +
scale_x_log10()
#Overlaying average marginal predicted probabilities on figure
#by Ref
smallDS_plot3 <- ggplot() +
geom_point(data = smallDS2, mapping = aes(
x = Sed_exposure_days,
y = Binary_small_necroses,
color = Ref)) +
labs(x = "Sediment exposure duration (days)",
y = "Predicted probability of small\nnecroses to sediment exposure",
color = "\nBinary Data by Study",
linetype = "Predicted\nProbability") +
scale_x_log10(limits=c(0.1,200), breaks=c(0.01,0.1,1,10,100,1000),
label=c("0.01","0.1","1","10","100","1000")) +
geom_ribbon(data = plotdat, aes(x = Sed_dur, y = PredictedProbability,
ymin = LowerQuantile, ymax = UpperQuantile),
alpha = 0.15) +
geom_line(data = plotdat, aes(x = Sed_dur, y = PredictedProbability,
linetype = "twodash")) +
geom_line(data = plotdat, aes(x = Sed_dur, y = MedianProbability,
linetype = "solid")) +
theme_bw() +
scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))
smallDS_plot3
## n
## 1 522
## boundary (singular) fit: see ?isSingular
## Data: largeDS2
## Models:
## modDSb6: Binary_large_necroses ~ Sed_level_stand_mg * Sed_exposure_days +
## modDSb6: (1 | Gsp)
## modDSb9: Binary_large_necroses ~ Sed_level_stand_mg * Sed_exposure_days +
## modDSb9: (1 | Ref)
## modDSb12: Binary_large_necroses ~ Sed_level_stand_mg * Sed_exposure_days +
## modDSb12: (1 | Gsp) + (1 | Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modDSb6 5 346.09 367.38 -168.04 336.09
## modDSb9 5 392.60 413.88 -191.30 382.60 0.000 0 1
## modDSb12 6 345.25 370.80 -166.62 333.25 49.346 1 2.146e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Data: largeDS2
## Models:
## modDSb6: Binary_large_necroses ~ Sed_level_stand_mg * Sed_exposure_days +
## modDSb6: (1 | Gsp)
## modDSb12: Binary_large_necroses ~ Sed_level_stand_mg * Sed_exposure_days +
## modDSb12: (1 | Gsp) + (1 | Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modDSb6 5 346.09 367.38 -168.04 336.09
## modDSb12 6 345.25 370.80 -166.62 333.25 2.8366 1 0.09214 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Data: largeDS2
## Models:
## modDSb4: Binary_large_necroses ~ Sed_level_stand_mg + (1 | Gsp)
## modDSb5: Binary_large_necroses ~ Sed_level_stand_mg + Sed_exposure_days +
## modDSb5: (1 | Gsp)
## modDSb6: Binary_large_necroses ~ Sed_level_stand_mg * Sed_exposure_days +
## modDSb6: (1 | Gsp)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modDSb4 3 350.54 363.31 -172.27 344.54
## modDSb5 4 345.35 362.38 -168.68 337.35 7.1909 1 0.007327 **
## modDSb6 5 346.09 367.38 -168.04 336.09 1.2629 1 0.261106
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Data: largeDS2
## Models:
## modDSb10: Binary_large_necroses ~ Sed_level_stand_mg + (1 | Gsp) + (1 |
## modDSb10: Ref)
## modDSb11: Binary_large_necroses ~ Sed_level_stand_mg + Sed_exposure_days +
## modDSb11: (1 | Gsp) + (1 | Ref)
## modDSb12: Binary_large_necroses ~ Sed_level_stand_mg * Sed_exposure_days +
## modDSb12: (1 | Gsp) + (1 | Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modDSb10 4 352.54 369.57 -172.27 344.54
## modDSb11 5 343.95 365.24 -166.98 333.95 10.589 1 0.001137 **
## modDSb12 6 345.25 370.80 -166.62 333.25 0.701 1 0.402432
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Data: largeDS2
## Models:
## modDSb5: Binary_large_necroses ~ Sed_level_stand_mg + Sed_exposure_days +
## modDSb5: (1 | Gsp)
## modDSb11: Binary_large_necroses ~ Sed_level_stand_mg + Sed_exposure_days +
## modDSb11: (1 | Gsp) + (1 | Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modDSb5 4 345.35 362.38 -168.68 337.35
## modDSb11 5 343.95 365.24 -166.98 333.95 3.3985 1 0.06526 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Binary_large_necroses ~ Sed_level_stand_mg + Sed_exposure_days +
## (1 | Gsp) + (1 | Ref)
## Data: largeDS2
##
## AIC BIC logLik deviance df.resid
## 344.0 365.2 -167.0 334.0 517
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9660 -0.3198 -0.0383 -0.0265 4.5830
##
## Random effects:
## Groups Name Variance Std.Dev.
## Gsp (Intercept) 28.3671 5.3261
## Ref (Intercept) 0.9538 0.9766
## Number of obs: 522, groups: Gsp, 74; Ref, 18
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -7.702214 2.000071 -3.851 0.000118 ***
## Sed_level_stand_mg 0.003963 0.001377 2.878 0.003998 **
## Sed_exposure_days 0.026231 0.009663 2.715 0.006635 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?
## Model is nearly unidentifiable: large eigenvalue ratio
## - Rescale variables?
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Binary_large_necroses ~ log10_conc + log10_dur + (1 | Gsp) +
## (1 | Ref)
## Data: largeDS2
##
## AIC BIC logLik deviance df.resid
## 322.9 344.2 -156.5 312.9 517
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.7377 -0.1869 -0.0426 -0.0201 4.1335
##
## Random effects:
## Groups Name Variance Std.Dev.
## Gsp (Intercept) 27.088 5.205
## Ref (Intercept) 1.146 1.071
## Number of obs: 522, groups: Gsp, 74; Ref, 18
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -14.0241 2.7141 -5.167 2.38e-07 ***
## log10_conc 2.5721 0.5767 4.460 8.20e-06 ***
## log10_dur 2.7644 0.8983 3.077 0.00209 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] 0.9099617
##
## p 0 1
## 0 434 34
## 1 13 41
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
## Area under the curve: 0.9494
## R2m R2c
## theoretical 0.09338887 0.9053851
## delta 0.09041040 0.8765094
## $Gsp
##
## $Ref
## Groups Name Std.Dev.
## Gsp (Intercept) 5.2046
## Ref (Intercept) 1.0707
## [1] 2.917504
## Est LL UL
## (Intercept) 8.117388e-07 3.973205e-09 1.658409e-04
## log10_conc 1.309276e+01 4.227962e+00 4.054444e+01
## log10_dur 1.586911e+01 2.728155e+00 9.230725e+01
For every 10-fold increase in exposure concentration of deposited sediment, the odds of any tissue mortality increase by 13.1 times (95% CI 4.2, 40.5; GLMM z = 4.460, p < 0.00001), while holding exposure duration constant and accounting for variability among studies and species. For every 10-fold change in exposure duration, the odds of any tissue mortality increase by 15.9 times (95% CI 2.7, 9.2; GLMM z = 3.077, p = 0.002), while holding exposure concentration constant and accounting for variability among studies and species.
#Overlaying glmm results on figure, but prediction line is jagged...
pred_modDSb11_log <- predict(modDSb11_log, newdata=largeDS2, type="response")
largeDS3 <- cbind(largeDS2, pred_modDSb11_log)
largeDS_plot <- ggplot(data = largeDS3) +
geom_point(mapping = aes(
x = Sed_level_stand_mg,
y = Binary_large_necroses,
color = Ref)) +
labs(x = expression("Sediment exposure concentration (mg/cm"^"2"*"/day)"),
y = "Large necroses due to sediment exposure?",
color = "Study") +
scale_x_log10(limits=c(1,1000),breaks=c(0.01,0.1,1,10,100,1000),
label=c("0.01","0.1","1","10","100","1000")) +
geom_line(aes(x = Sed_level_stand_mg, y = pred_modDSb11_log), inherit.aes=FALSE) +
theme_bw()
largeDS_plot
That plot is confusing to interpret. Let’s see if I can plot the average marginal probability, i.e., the average change in probability of the outcome across the range of the predictor of interest. This is described in some detail at the following, useful website: https://stats.idre.ucla.edu/r/dae/mixed-effects-logistic-regression/
jvalues <- with(largeDS2, seq(from = min(log10_conc), to = max(log10_conc), length.out = 100))
# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
largeDS2$log10_conc <- j
predict(modDSb11_log, newdata = largeDS2, type = "response")
})
# average marginal predicted probability across a few different sediment levels
#sapply(pp[c(1, 20, 40, 60, 80, 100)], mean)
##
# get the means with lower and upper quartiles
plotdat <- t(sapply(pp, function(x) {
c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))
# add in log10_conc values and convert to data frame
plotdat <- as.data.frame(cbind(plotdat, jvalues))
plotdat <- plotdat %>% mutate(Sed_level_linear=10^jvalues)
# better names and show the first few rows
colnames(plotdat) <- c("PredictedProbability", "MedianProbability",
"LowerQuantile", "UpperQuantile", "log10_conc", "Sed_conc")
#head(plotdat)
# plot average marginal predicted probabilities
ggplot(plotdat, aes(x = log10_conc)) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
ggplot(plotdat, aes(x = log10_conc)) +
geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
ggplot(plotdat, aes(x = Sed_conc, y = PredictedProbability)) +
geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1)) +
scale_x_log10()
#Overlaying average marginal predicted probabilities on figure
#by Ref
largeDS_plot2 <- ggplot() +
geom_point(data = largeDS2, mapping = aes(
x = Sed_level_stand_mg,
y = Binary_large_necroses)) +
labs(x = expression("Sediment exposure concentration (mg/cm"^"2"*"/day)"),
y = "Predicted probability of large\nnecroses due to sediment exposure",
linetype = "Predicted\nProbability") +
scale_x_log10(limits=c(1,1000), breaks=c(0.01,0.1,1,10,100,1000,loaelDS2),
label=c("0.01","0.1","1","10","100","1000",round(loaelDS2,digits=1))) +
geom_ribbon(data = plotdat, aes(x = Sed_conc, y = PredictedProbability,
ymin = LowerQuantile, ymax = UpperQuantile),
alpha = 0.15) +
geom_line(data = plotdat, aes(x = Sed_conc, y = PredictedProbability,
linetype = "twodash")) +
geom_line(data = plotdat, aes(x = Sed_conc, y = MedianProbability,
linetype = "solid")) +
geom_vline(xintercept=loaelDS2, linetype="dashed", color = "red") +
theme_bw() +
scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))
largeDS_plot2
jvalues <- with(largeDS2, seq(from = min(log10_dur), to = max(log10_dur), length.out = 100))
# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
largeDS2$log10_dur <- j
predict(modDSb11_log, newdata = largeDS2, type = "response")
})
# average marginal predicted probability across a few different sediment levels
#sapply(pp[c(1, 20, 40, 60, 80, 100)], mean)
##
# get the means with lower and upper quartiles
plotdat <- t(sapply(pp, function(x) {
c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))
# add in log10_dur values and convert to data frame
plotdat <- as.data.frame(cbind(plotdat, jvalues))
plotdat <- plotdat %>% mutate(Sed_level_linear=10^jvalues)
# better names and show the first few rows
colnames(plotdat) <- c("PredictedProbability", "MedianProbability",
"LowerQuantile", "UpperQuantile", "log10_dur", "Sed_dur")
#head(plotdat)
# plot average marginal predicted probabilities
ggplot(plotdat, aes(x = log10_dur)) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
ggplot(plotdat, aes(x = log10_dur)) +
geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
ggplot(plotdat, aes(x = Sed_dur, y = PredictedProbability)) +
geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1)) +
scale_x_log10()
#Overlaying average marginal predicted probabilities on figure
#by Ref
largeDS_plot2 <- ggplot() +
geom_point(data = largeDS2, mapping = aes(
x = Sed_exposure_days,
y = Binary_large_necroses)) +
labs(x = expression("Sediment exposure duration (days)"),
y = "Predicted probability of large\nnecroses due to sediment exposure",
linetype = "Predicted\nProbability") +
scale_x_log10(limits=c(0.1,max(largeDS2$Sed_exposure_days)), breaks=c(0.01,0.1,1,10,100,1000),
label=c("0.01","0.1","1","10","100","1000")) +
geom_ribbon(data = plotdat, aes(x = Sed_dur, y = PredictedProbability,
ymin = LowerQuantile, ymax = UpperQuantile),
alpha = 0.15) +
geom_line(data = plotdat, aes(x = Sed_dur, y = PredictedProbability,
linetype = "twodash")) +
geom_line(data = plotdat, aes(x = Sed_dur, y = MedianProbability,
linetype = "solid")) +
theme_bw() +
scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))
largeDS_plot2
## n
## 1 509
## Data: totalDS2
## Models:
## modDSc6: Binary_death_total ~ Sed_level_stand_mg * Sed_exposure_days +
## modDSc6: (1 | Gsp)
## modDSc9: Binary_death_total ~ Sed_level_stand_mg * Sed_exposure_days +
## modDSc9: (1 | Ref)
## modDSc12: Binary_death_total ~ Sed_level_stand_mg * Sed_exposure_days +
## modDSc12: (1 | Gsp) + (1 | Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modDSc6 5 306.66 327.82 -148.33 296.66
## modDSc9 5 392.75 413.91 -191.37 382.75 0.000 0 1
## modDSc12 6 295.74 321.13 -141.87 283.74 99.009 1 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Data: totalDS2
## Models:
## modDSc10: Binary_death_total ~ Sed_level_stand_mg + (1 | Gsp) + (1 | Ref)
## modDSc11: Binary_death_total ~ Sed_level_stand_mg + Sed_exposure_days +
## modDSc11: (1 | Gsp) + (1 | Ref)
## modDSc12: Binary_death_total ~ Sed_level_stand_mg * Sed_exposure_days +
## modDSc12: (1 | Gsp) + (1 | Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modDSc10 4 308.19 325.12 -150.09 300.19
## modDSc11 5 295.52 316.69 -142.76 285.52 14.664 1 0.0001285 ***
## modDSc12 6 295.74 321.13 -141.87 283.74 1.784 1 0.1816564
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Data: totalDS2
## Models:
## modDSc4: Binary_death_total ~ Sed_level_stand_mg + (1 | Gsp)
## modDSc7: Binary_death_total ~ Sed_level_stand_mg + (1 | Ref)
## modDSc10: Binary_death_total ~ Sed_level_stand_mg + (1 | Gsp) + (1 | Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modDSc4 3 319.21 331.90 -156.60 313.21
## modDSc7 3 396.32 409.02 -195.16 390.32 0.000 0 1
## modDSc10 4 308.19 325.12 -150.09 300.19 90.131 1 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Binary_death_total ~ log10_conc + log10_dur + (1 | Gsp) + (1 |
## Ref)
## Data: totalDS2
##
## AIC BIC logLik deviance df.resid
## 255.0 276.2 -122.5 245.0 504
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2021 -0.0137 -0.0021 0.0000 7.6530
##
## Random effects:
## Groups Name Variance Std.Dev.
## Gsp (Intercept) 280.27 16.741
## Ref (Intercept) 42.47 6.517
## Number of obs: 509, groups: Gsp, 75; Ref, 23
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -38.909 6.755 -5.760 8.42e-09 ***
## log10_conc 6.103 1.500 4.069 4.73e-05 ***
## log10_dur 7.855 1.711 4.590 4.44e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] 0.9685658
##
## p 0 1
## 0 419 9
## 1 7 74
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
## Area under the curve: 0.9905
## R2m R2c
## theoretical 0.0683368385 0.99059879
## delta 0.0009323897 0.01351576
## $Gsp
##
## $Ref
## Groups Name Std.Dev.
## Gsp (Intercept) 16.741
## Ref (Intercept) 6.517
## [1] 676.5725
## Est LL UL
## (Intercept) 1.264562e-17 2.247393e-23 7.115428e-12
## log10_conc 4.471818e+02 2.363854e+01 8.459556e+03
## log10_dur 2.577456e+03 9.003714e+01 7.378375e+04
For every 10-fold increase in exposure concentration of deposited sediment, the odds of total colony mortality increase by 447 times (95% CI 24, 8459; GLMM z = 4.069, p < 0.0001), after accounting for exposure duration and variability among studies and species. For every 10-fold change in exposure duration, the odds of total colony mortality increase by 2577 times (95% CI 90, 73,784; GLMM z = 4.590, p = 0.00001), after accounting for exposure concentration and variability among studies and species.
#Overlaying glmm results on figure, but prediction line is jagged...
pred_modDSc11_log <- predict(modDSc11_log, newdata=totalDS2, type="response")
totalDS3 <- cbind(totalDS2, pred_modDSc11_log)
totalDS_plot <- ggplot(data = totalDS3) +
geom_point(mapping = aes(
x = Sed_level_stand_mg,
y = Binary_death_total,
color = Ref)) +
labs(x = expression("Sediment exposure concentration (mg/cm"^"2"*"/day)"),
y = "Total mortality due to sediment exposure?",
color = "Study") +
scale_x_log10(limits=c(0.01,1000),breaks=c(0.01,0.1,1,10,100,1000),
label=c("0.01","0.1","1","10","100","1000")) +
geom_line(aes(x = Sed_level_stand_mg, y = pred_modDSc11_log), inherit.aes=FALSE) +
theme_bw()
totalDS_plot
That plot is confusing to interpret. Let’s see if I can plot the average marginal probability, i.e., the average change in probability of the outcome across the range of the predictor of interest. This is described in some detail at the following, useful website: https://stats.idre.ucla.edu/r/dae/mixed-effects-logistic-regression/
jvalues <- with(totalDS2, seq(from = min(log10_conc), to = max(log10_conc), length.out = 100))
# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
totalDS2$log10_conc <- j
predict(modDSc11_log, newdata = totalDS2, type = "response")
})
# average marginal predicted probability across a few different sediment levels
#sapply(pp[c(1, 20, 40, 60, 80, 100)], mean)
##
# get the means with lower and upper quartiles
plotdat <- t(sapply(pp, function(x) {
c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))
# add in log10_conc values and convert to data frame
plotdat <- as.data.frame(cbind(plotdat, jvalues))
plotdat <- plotdat %>% mutate(Sed_level_linear=10^jvalues)
# better names and show the first few rows
colnames(plotdat) <- c("PredictedProbability", "MedianProbability",
"LowerQuantile", "UpperQuantile", "log10_conc", "Sed_conc")
#head(plotdat)
# plot average marginal predicted probabilities
ggplot(plotdat, aes(x = log10_conc)) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
ggplot(plotdat, aes(x = log10_conc)) +
geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
ggplot(plotdat, aes(x = Sed_conc, y = PredictedProbability)) +
geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1)) +
scale_x_log10()
#Overlaying average marginal predicted probabilities on figure
#by Ref
totalDS_plot2 <- ggplot() +
geom_point(data = totalDS2, mapping = aes(
x = Sed_level_stand_mg,
y = Binary_death_total,
color = Ref)) +
labs(x = expression("Sediment exposure concentration (mg/cm"^"2"*"/day)"),
y = "Predicted probability of total colony\nmortality due to sediment exposure",
color = "Study",
linetype = "Predicted Probability") +
scale_x_log10(limits=c(0.01,1000), breaks=c(0.01,0.1,1,10,100,1000,loaelDS3),
label=c("0.01","0.1","1","10","100","1000",round(loaelDS3,digits=1))) +
geom_ribbon(data = plotdat, aes(x = Sed_conc, y = PredictedProbability,
ymin = LowerQuantile, ymax = UpperQuantile),
alpha = 0.15) +
geom_line(data = plotdat, aes(x = Sed_conc, y = PredictedProbability,
linetype = "twodash")) +
geom_line(data = plotdat, aes(x = Sed_conc, y = MedianProbability,
linetype = "solid")) +
geom_vline(xintercept=loaelDS3, linetype="dashed", color = "red") +
theme_bw() +
scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))
totalDS_plot2
jvalues <- with(totalDS2, seq(from = min(log10_dur), to = max(log10_dur), length.out = 100))
# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
totalDS2$log10_dur <- j
predict(modDSc11_log, newdata = totalDS2, type = "response")
})
# get the means with lower and upper quartiles
plotdat <- t(sapply(pp, function(x) {
c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))
# add in log10_dur values and convert to data frame
plotdat <- as.data.frame(cbind(plotdat, jvalues))
plotdat <- plotdat %>% mutate(Sed_level_linear=10^jvalues)
# better names and show the first few rows
colnames(plotdat) <- c("PredictedProbability", "MedianProbability",
"LowerQuantile", "UpperQuantile", "log10_dur", "Sed_dur")
#head(plotdat)
# plot average marginal predicted probabilities
ggplot(plotdat, aes(x = log10_dur)) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
ggplot(plotdat, aes(x = log10_dur)) +
geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
ggplot(plotdat, aes(x = Sed_dur, y = PredictedProbability)) +
geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1)) +
scale_x_log10()
#Overlaying average marginal predicted probabilities on figure
#by Ref
totalDS_plot2 <- ggplot() +
geom_point(data = totalDS2, mapping = aes(
x = Sed_exposure_days,
y = Binary_death_total)) +
labs(x = expression("Sediment exposure duration (days)"),
y = "Predicted probability of total colony\nmortality due to sediment exposure",
linetype = "Predicted\nProbability") +
scale_x_log10(limits=c(0.1,max(totalDS2$Sed_exposure_days)), breaks=c(0.01,0.1,1,10,100,1000),
label=c("0.01","0.1","1","10","100","1000")) +
geom_ribbon(data = plotdat, aes(x = Sed_dur, y = PredictedProbability,
ymin = LowerQuantile, ymax = UpperQuantile),
alpha = 0.15) +
geom_line(data = plotdat, aes(x = Sed_dur, y = PredictedProbability,
linetype = "twodash")) +
geom_line(data = plotdat, aes(x = Sed_dur, y = MedianProbability,
linetype = "solid")) +
theme_bw() +
scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))
totalDS_plot2
## n
## 1 147
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## Data: smallSS2
## Models:
## modSS6: Binary_small_necroses ~ Sed_level_stand_mg * Sed_exposure_days +
## modSS6: (1 | Gsp)
## modSS9: Binary_small_necroses ~ Sed_level_stand_mg * Sed_exposure_days +
## modSS9: (1 | Ref)
## modSS12: Binary_small_necroses ~ Sed_level_stand_mg * Sed_exposure_days +
## modSS12: (1 | Gsp) + (1 | Ref)
## modSS15: Binary_small_necroses ~ Sed_level_stand_mg * Sed_exposure_days +
## modSS15: (1 | Ref_name/Ref)
## modSS18: Binary_small_necroses ~ Sed_level_stand_mg * Sed_exposure_days +
## modSS18: (1 | Gsp) + (1 | Ref_name)
## modSS21: Binary_small_necroses ~ Sed_level_stand_mg * Sed_exposure_days +
## modSS21: (1 | Gsp) + (1 | Ref_name/Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modSS6 5 108.19 123.14 -49.094 98.189
## modSS9 5 126.29 141.24 -58.145 116.290 0.000 0 1
## modSS12 6 110.19 128.13 -49.094 98.189 18.102 1 2.094e-05 ***
## modSS15 6 128.29 146.23 -58.145 116.290 0.000 0 1
## modSS18 6 110.19 128.13 -49.094 98.189 18.102 0 < 2.2e-16 ***
## modSS21 7 112.19 133.12 -49.094 98.189 0.000 1 1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Data: smallSS2
## Models:
## modSS6: Binary_small_necroses ~ Sed_level_stand_mg * Sed_exposure_days +
## modSS6: (1 | Gsp)
## modSS12: Binary_small_necroses ~ Sed_level_stand_mg * Sed_exposure_days +
## modSS12: (1 | Gsp) + (1 | Ref)
## modSS18: Binary_small_necroses ~ Sed_level_stand_mg * Sed_exposure_days +
## modSS18: (1 | Gsp) + (1 | Ref_name)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modSS6 5 108.19 123.14 -49.094 98.189
## modSS12 6 110.19 128.13 -49.094 98.189 0 1 1
## modSS18 6 110.19 128.13 -49.094 98.189 0 0 1
## Data: smallSS2
## Models:
## modSS4: Binary_small_necroses ~ Sed_level_stand_mg + (1 | Gsp)
## modSS5: Binary_small_necroses ~ Sed_level_stand_mg + Sed_exposure_days +
## modSS5: (1 | Gsp)
## modSS6: Binary_small_necroses ~ Sed_level_stand_mg * Sed_exposure_days +
## modSS6: (1 | Gsp)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modSS4 3 135.52 144.49 -64.759 129.518
## modSS5 4 106.75 118.71 -49.375 98.750 30.7680 1 2.908e-08 ***
## modSS6 5 108.19 123.14 -49.094 98.189 0.5607 1 0.454
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Data: smallSS2
## Models:
## modSS10: Binary_small_necroses ~ Sed_level_stand_mg + (1 | Gsp) + (1 |
## modSS10: Ref)
## modSS11: Binary_small_necroses ~ Sed_level_stand_mg + Sed_exposure_days +
## modSS11: (1 | Gsp) + (1 | Ref)
## modSS12: Binary_small_necroses ~ Sed_level_stand_mg * Sed_exposure_days +
## modSS12: (1 | Gsp) + (1 | Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modSS10 4 137.45 149.41 -64.724 129.447
## modSS11 5 108.75 123.70 -49.375 98.750 30.6976 1 3.015e-08 ***
## modSS12 6 110.19 128.13 -49.094 98.189 0.5607 1 0.454
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Data: smallSS2
## Models:
## modSS16: Binary_small_necroses ~ Sed_level_stand_mg + (1 | Gsp) + (1 |
## modSS16: Ref_name)
## modSS17: Binary_small_necroses ~ Sed_level_stand_mg + Sed_exposure_days +
## modSS17: (1 | Gsp) + (1 | Ref_name)
## modSS18: Binary_small_necroses ~ Sed_level_stand_mg * Sed_exposure_days +
## modSS18: (1 | Gsp) + (1 | Ref_name)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modSS16 4 137.45 149.41 -64.724 129.447
## modSS17 5 108.75 123.70 -49.375 98.750 30.6976 1 3.015e-08 ***
## modSS18 6 110.19 128.13 -49.094 98.189 0.5607 1 0.454
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Data: smallSS2
## Models:
## modSS5: Binary_small_necroses ~ Sed_level_stand_mg + Sed_exposure_days +
## modSS5: (1 | Gsp)
## modSS11: Binary_small_necroses ~ Sed_level_stand_mg + Sed_exposure_days +
## modSS11: (1 | Gsp) + (1 | Ref)
## modSS17: Binary_small_necroses ~ Sed_level_stand_mg + Sed_exposure_days +
## modSS17: (1 | Gsp) + (1 | Ref_name)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modSS5 4 106.75 118.71 -49.375 98.75
## modSS11 5 108.75 123.70 -49.375 98.75 0 1 1
## modSS17 5 108.75 123.70 -49.375 98.75 0 0 1
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Binary_small_necroses ~ Sed_level_stand_mg + Sed_exposure_days +
## (1 | Gsp)
## Data: smallSS2
##
## AIC BIC logLik deviance df.resid
## 106.7 118.7 -49.4 98.7 143
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.9893 -0.3421 -0.1326 -0.0291 3.8738
##
## Random effects:
## Groups Name Variance Std.Dev.
## Gsp (Intercept) 3.44 1.855
## Number of obs: 147, groups: Gsp, 8
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -8.06121 1.75339 -4.598 4.28e-06 ***
## Sed_level_stand_mg 0.03325 0.00998 3.331 0.000864 ***
## Sed_exposure_days 0.07866 0.01806 4.355 1.33e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] 0.8707483
##
## p 0 1
## 0 108 11
## 1 8 20
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
## Area under the curve: 0.9189
## R2m R2c
## theoretical 0.3536090 0.6839914
## delta 0.2701203 0.5224978
## $Gsp
## Groups Name Std.Dev.
## Gsp (Intercept) 1.8546
## [1] 6.389086
## Est LL UL
## (Intercept) 0.003743975 0.0003457821 0.04053811
## Sed_level_stand_mg 1.023311817 1.0095310571 1.03728069
## Sed_exposure_days 1.056033682 1.0304395187 1.08226356
## Est LL UL
## (Intercept) 8.685337e-09 3.177978e-12 2.373682e-05
## Sed_level_stand_mg 1.079558e+00 1.032013e+00 1.129292e+00
## Sed_exposure_days 1.198549e+00 1.104739e+00 1.300325e+00
For every 10-fold increase in exposure concentration of suspended sediment, the odds of developing small tissue necroses increase by 1.08 times (95% CI 1.03, 1.13; GLMM z = 3.331, p = 0.0008), while holding exposure duration constant and after accounting for variability among species. For every 10-fold increase in exposure duration, the odds of small tissue necroses increase by 1.20 times (95% CI 1.10, 1.30; GLMM z = 4.355, p < 0.0001), while holding exposure concentration constant and after accounting for variability among species.
#Overlaying glmm results on figure, but prediction line is jagged...
pred_modSS5 <- predict(modSS5, newdata=smallSS2, type="response")
smallSS3 <- cbind(smallSS2, pred_modSS5)
smallSS_plot <- ggplot(data = smallSS3) +
geom_point(mapping = aes(
x = Sed_level_stand_mg,
y = Binary_small_necroses,
color = Gsp)) +
labs(x = "Sediment exposure concentration (mg/L)",
y = "Small necroses due to sediment exposure?",
color = "Species") +
scale_x_log10(limits=c(1,max(smallSS2$Sed_level_stand_mg)),breaks=c(0.01,0.1,1,10,100,1000),
label=c("0.01","0.1","1","10","100","1000")) +
geom_line(aes(x = Sed_level_stand_mg, y = pred_modSS5), inherit.aes=FALSE) +
theme_bw()
smallSS_plot
That plot is confusing to interpret. Let’s see if I can plot the average marginal probability, i.e., the average change in probability of the outcome across the range of the predictor of interest. This is described in some detail at the following, useful website: https://stats.idre.ucla.edu/r/dae/mixed-effects-logistic-regression/
jvalues <- with(smallSS2, seq(from = min(Sed_level_stand_mg), to = max(Sed_level_stand_mg), length.out = 100))
# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
smallSS2$Sed_level_stand_mg <- j
predict(modSS5, newdata = smallSS2, type = "response")
})
# get the means with lower and upper quartiles
plotdat <- t(sapply(pp, function(x) {
c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))
# add in Sed_level_stand_mg values and convert to data frame
plotdat <- as.data.frame(cbind(plotdat, jvalues))
# better names and show the first few rows
colnames(plotdat) <- c("PredictedProbability", "MedianProbability",
"LowerQuantile", "UpperQuantile", "Sed_conc")
#head(plotdat)
# plot average marginal predicted probabilities
ggplot(plotdat, aes(x = Sed_conc)) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
ggplot(plotdat, aes(x = Sed_conc)) +
geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
#Overlaying average marginal predicted probabilities on figure
#by Ref
smallSS_plot2 <- ggplot() +
geom_point(data = smallSS2, mapping = aes(
x = Sed_level_stand_mg,
y = Binary_small_necroses)) +
labs(x = "Sediment exposure concentration (mg/L)",
y = "Predicted probability of small\nnecroses due to sediment exposure",
linetype = "Predicted\nProbability") +
scale_x_continuous(limits=c(0,125), breaks=c(0,25,50,75,100,125,loaelSS1),
labels=c("0","25","50","75","100","125",round(loaelSS1,digits=1))) +
geom_ribbon(data = plotdat, aes(x = Sed_conc, y = PredictedProbability,
ymin = LowerQuantile, ymax = UpperQuantile),
alpha = 0.15) +
geom_line(data = plotdat, aes(x = Sed_conc, y = PredictedProbability,
linetype = "twodash")) +
geom_line(data = plotdat, aes(x = Sed_conc, y = MedianProbability,
linetype = "solid")) +
geom_vline(xintercept=loaelSS1, linetype="dashed", color = "red") +
theme_bw() +
scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))
smallSS_plot2
jvalues <- with(smallSS2, seq(from = min(Sed_exposure_days), to = max(Sed_exposure_days), length.out = 100))
# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
smallSS2$Sed_exposure_days <- j
predict(modSS5, newdata = smallSS2, type = "response")
})
# get the means with lower and upper quartiles
plotdat <- t(sapply(pp, function(x) {
c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))
# add in Sed_exposure_days values and convert to data frame
plotdat <- as.data.frame(cbind(plotdat, jvalues))
# better names and show the first few rows
colnames(plotdat) <- c("PredictedProbability", "MedianProbability",
"LowerQuantile", "UpperQuantile", "Sed_dur")
#head(plotdat)
# plot average marginal predicted probabilities
ggplot(plotdat, aes(x = Sed_dur)) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
ggplot(plotdat, aes(x = Sed_dur)) +
geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
#Overlaying average marginal predicted probabilities on figure
#by Ref
smallSS_plot3 <- ggplot() +
geom_point(data = smallSS2, mapping = aes(
x = Sed_exposure_days,
y = Binary_small_necroses,
color = Gsp)) +
labs(x = "Sediment exposure duration (days)",
y = "Predicted probability of small\nnecroses to sediment exposure",
color = "Species",
linetype = "Predicted\nProbability") +
scale_x_continuous(limits=c(0,max(smallSS2$Sed_exposure_days))) +
geom_ribbon(data = plotdat, aes(x = Sed_dur, y = PredictedProbability,
ymin = LowerQuantile, ymax = UpperQuantile),
alpha = 0.15) +
geom_line(data = plotdat, aes(x = Sed_dur, y = PredictedProbability,
linetype = "twodash")) +
geom_line(data = plotdat, aes(x = Sed_dur, y = MedianProbability,
linetype = "solid")) +
theme_bw() +
scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))
smallSS_plot3
## n
## 1 147
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## df AIC
## modSSb1 2 52.19214
## modSSb2 3 11.54518
## modSSb3 4 13.54518
## Data: largeSS2
## Models:
## modSSb7: Binary_large_necroses ~ Sed_level_stand_mg + (1 | Ref)
## modSSb8: Binary_large_necroses ~ Sed_level_stand_mg + Sed_exposure_days +
## modSSb8: (1 | Ref)
## modSSb9: Binary_large_necroses ~ Sed_level_stand_mg * Sed_exposure_days +
## modSSb9: (1 | Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modSSb7 3 31.918 40.889 -12.9588 25.9176
## modSSb8 4 13.545 25.507 -2.7726 5.5452 20.372 1 6.374e-06 ***
## modSSb9 5 15.545 30.497 -2.7726 5.5452 0.000 1 1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Data: largeSS2
## Models:
## modSSb13: Binary_large_necroses ~ Sed_level_stand_mg + (1 | Ref_name/Ref)
## modSSb14: Binary_large_necroses ~ Sed_level_stand_mg + Sed_exposure_days +
## modSSb14: (1 | Ref_name/Ref)
## modSSb15: Binary_large_necroses ~ Sed_level_stand_mg * Sed_exposure_days +
## modSSb15: (1 | Ref_name/Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modSSb13 4 33.918 45.879 -12.9588 25.9176
## modSSb14 5 15.545 30.497 -2.7726 5.5452 20.372 1 6.374e-06 ***
## modSSb15 6 17.545 35.488 -2.7726 5.5452 0.000 1 1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Data: largeSS2
## Models:
## modSSb4: Binary_large_necroses ~ Sed_level_stand_mg + (1 | Gsp)
## modSSb8: Binary_large_necroses ~ Sed_level_stand_mg + Sed_exposure_days +
## modSSb8: (1 | Ref)
## modSSb10: Binary_large_necroses ~ Sed_level_stand_mg + (1 | Gsp) + (1 |
## modSSb10: Ref)
## modSSb16: Binary_large_necroses ~ Sed_level_stand_mg + (1 | Gsp) + (1 |
## modSSb16: Ref_name)
## modSSb14: Binary_large_necroses ~ Sed_level_stand_mg + Sed_exposure_days +
## modSSb14: (1 | Ref_name/Ref)
## modSSb19: Binary_large_necroses ~ Sed_level_stand_mg + (1 | Gsp) + (1 |
## modSSb19: Ref_name/Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modSSb4 3 36.784 45.756 -15.3921 30.7842
## modSSb8 4 13.545 25.507 -2.7726 5.5452 25.239 1 5.065e-07 ***
## modSSb10 4 33.803 45.765 -12.9014 25.8028 0.000 0 1
## modSSb16 4 33.803 45.765 -12.9014 25.8028 0.000 0 1
## modSSb14 5 15.545 30.497 -2.7726 5.5452 20.258 1 6.768e-06 ***
## modSSb19 5 35.803 50.755 -12.9014 25.8028 0.000 0 1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Data: largeSS2
## Models:
## modSSb8: Binary_large_necroses ~ Sed_level_stand_mg + Sed_exposure_days +
## modSSb8: (1 | Ref)
## modSSb14: Binary_large_necroses ~ Sed_level_stand_mg + Sed_exposure_days +
## modSSb14: (1 | Ref_name/Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modSSb8 4 13.545 25.507 -2.7726 5.5452
## modSSb14 5 15.545 30.497 -2.7726 5.5452 0 1 1
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Binary_large_necroses ~ Sed_level_stand_mg + Sed_exposure_days +
## (1 | Ref)
## Data: largeSS2
##
## AIC BIC logLik deviance df.resid
## 13.5 25.5 -2.8 5.5 143
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1 0 0 0 1
##
## Random effects:
## Groups Name Variance Std.Dev.
## Ref (Intercept) 0 0
## Number of obs: 147, groups: Ref, 4
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -779.798 968.212 -0.805 0.4206
## Sed_level_stand_mg 2.026 31.654 0.064 0.9490
## Sed_exposure_days 8.582 3.984 2.154 0.0312 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Binary_large_necroses ~ log10_conc + log10_dur + (1 | Ref)
## Data: largeSS2
##
## AIC BIC logLik deviance df.resid
## 13.5 25.5 -2.8 5.5 143
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1 0 0 0 1
##
## Random effects:
## Groups Name Variance Std.Dev.
## Ref (Intercept) 3.193e-17 5.65e-09
## Number of obs: 147, groups: Ref, 4
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1347.0 49256351.2 0 1
## log10_conc 85.8 16165394.5 0 1
## log10_dur 634.7 15817710.9 0 1
## convergence code: 0
## boundary (singular) fit: see ?isSingular
##
## Call:
## glm(formula = Binary_large_necroses ~ Sed_level_stand_mg + Sed_exposure_days,
## family = binomial(link = "logit"), data = largeSS2)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.177 0.000 0.000 0.000 1.177
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -390.084 100469.392 -0.004 0.997
## Sed_level_stand_mg 1.077 314.538 0.003 0.997
## Sed_exposure_days 4.271 1100.945 0.004 0.997
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 50.1358 on 146 degrees of freedom
## Residual deviance: 5.5452 on 144 degrees of freedom
## AIC: 11.545
##
## Number of Fisher Scoring iterations: 25
## [1] 0.9863946
##
## p 0 1
## 0 141 2
## 1 0 4
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
## Area under the curve: 0.9976
## R2m R2c
## theoretical 0.9996563 0.9996563
## delta 0.9973381 0.9973381
There is no significant relationship between sediment exposure and the odds of developing large tissue necroses (concentration: GLMM z = 0.003, p = 0.997; duration: GLMM z = 0.004, p = 0.997).
#Overlaying glmm results on figure, but prediction line is jagged...
pred_modSSb2 <- predict(modSSb2, newdata=largeSS2, type="response")
largeSS3 <- cbind(largeSS2, pred_modSSb2)
largeSS_plot <- ggplot(data = largeSS3) +
geom_point(mapping = aes(
x = Sed_level_stand_mg,
y = Binary_large_necroses,
color = Ref)) +
labs(x = "Sediment exposure concentration (mg/L)",
y = "Large necroses due to sediment exposure?",
color = "Study") +
scale_x_continuous(limits=c(0,max(largeSS2$Sed_level_stand_mg))) +
geom_line(aes(x = Sed_level_stand_mg, y = pred_modSSb2), inherit.aes=FALSE) +
theme_bw()
largeSS_plot
That plot is confusing to interpret. Let’s see if I can plot the average marginal probability, i.e., the average change in probability of the outcome across the range of the predictor of interest. This is described in some detail at the following, useful website: https://stats.idre.ucla.edu/r/dae/mixed-effects-logistic-regression/
jvalues <- with(largeSS2, seq(from = min(Sed_level_stand_mg), to = max(Sed_level_stand_mg), length.out = 100))
# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
largeSS2$Sed_level_stand_mg <- j
predict(modSSb2, newdata = largeSS2, type = "response")
})
# get the means with lower and upper quartiles
plotdat <- t(sapply(pp, function(x) {
c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))
# add in Sed_level_stand_mg values and convert to data frame
plotdat <- as.data.frame(cbind(plotdat, jvalues))
# better names and show the first few rows
colnames(plotdat) <- c("PredictedProbability", "MedianProbability",
"LowerQuantile", "UpperQuantile", "Sed_conc")
#head(plotdat)
# plot average marginal predicted probabilities
ggplot(plotdat, aes(x = Sed_conc)) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
ggplot(plotdat, aes(x = Sed_conc)) +
geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
#Overlaying average marginal predicted probabilities on figure
#by Ref
largeSS_plot2 <- ggplot() +
geom_point(data = largeSS2, mapping = aes(
x = Sed_level_stand_mg,
y = Binary_large_necroses,
color = Ref)) +
labs(x = "Sediment exposure concentration (mg/L)",
y = "Predicted probability of large\nnecroses due to sediment exposure",
color = "Binary data\nby Study",
linetype = "Predicted\nProbability") +
scale_x_continuous(limits=c(0,max(largeSS2$Sed_level_stand_mg)),
breaks = c(0,25,50,75,100,125,loaelSS2),
labels = c("0","25","50","75","100","125",round(loaelSS2,digits=1))) +
geom_ribbon(data = plotdat, aes(x = Sed_conc, y = PredictedProbability,
ymin = LowerQuantile, ymax = UpperQuantile),
alpha = 0.15) +
geom_line(data = plotdat, aes(x = Sed_conc, y = PredictedProbability,
linetype = "twodash")) +
geom_line(data = plotdat, aes(x = Sed_conc, y = MedianProbability,
linetype = "solid")) +
geom_vline(xintercept=loaelSS2, linetype="dashed", color = "red") +
theme_bw() +
scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))
largeSS_plot2
jvalues <- with(largeSS2, seq(from = min(Sed_exposure_days), to = max(Sed_exposure_days), length.out = 100))
# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
largeSS2$Sed_exposure_days <- j
predict(modSSb2, newdata = largeSS2, type = "response")
})
# get the means with lower and upper quartiles
plotdat <- t(sapply(pp, function(x) {
c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))
# add in Sed_exposure_days values and convert to data frame
plotdat <- as.data.frame(cbind(plotdat, jvalues))
# better names and show the first few rows
colnames(plotdat) <- c("PredictedProbability", "MedianProbability",
"LowerQuantile", "UpperQuantile", "Sed_dur")
#head(plotdat)
# plot average marginal predicted probabilities
ggplot(plotdat, aes(x = Sed_dur)) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
ggplot(plotdat, aes(x = Sed_dur)) +
geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
#Overlaying average marginal predicted probabilities on figure
#by Ref
largeSS_plot2 <- ggplot() +
geom_point(data = largeSS2, mapping = aes(
x = Sed_exposure_days,
y = Binary_large_necroses)) +
labs(x = expression("Sediment exposure duration (days)"),
y = "Predicted probability of large\nnecroses due to sediment exposure",
linetype = "Predicted\nProbability") +
geom_ribbon(data = plotdat, aes(x = Sed_dur, y = PredictedProbability,
ymin = LowerQuantile, ymax = UpperQuantile),
alpha = 0.15) +
geom_line(data = plotdat, aes(x = Sed_dur, y = PredictedProbability,
linetype = "twodash")) +
geom_line(data = plotdat, aes(x = Sed_dur, y = MedianProbability,
linetype = "solid")) +
theme_bw() +
scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))
largeSS_plot2
## n
## 1 176
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## df AIC
## modSSc1 2 56.11582
## modSSc2 3 11.54518
## modSSc3 4 13.54518
## Data: totalSS2
## Models:
## modSSc7: Binary_death_total ~ Sed_level_stand_mg + (1 | Ref)
## modSSc8: Binary_death_total ~ Sed_level_stand_mg + Sed_exposure_days +
## modSSc8: (1 | Ref)
## modSSc9: Binary_death_total ~ Sed_level_stand_mg * Sed_exposure_days +
## modSSc9: (1 | Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modSSc7 3 32.442 41.953 -13.2208 26.4416
## modSSc8 4 13.545 26.227 -2.7726 5.5452 20.896 1 4.848e-06 ***
## modSSc9 5 15.545 31.398 -2.7726 5.5452 0.000 1 1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Data: totalSS2
## Models:
## modSSc13: Binary_death_total ~ Sed_level_stand_mg + (1 | Ref_name/Ref)
## modSSc15: Binary_death_total ~ Sed_level_stand_mg * Sed_exposure_days +
## modSSc15: (1 | Ref_name/Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modSSc13 4 34.441 47.122 -13.2203 26.4406
## modSSc15 6 17.545 36.568 -2.7726 5.5452 20.895 2 2.902e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Data: totalSS2
## Models:
## modSSc4: Binary_death_total ~ Sed_level_stand_mg + (1 | Gsp)
## modSSc8: Binary_death_total ~ Sed_level_stand_mg + Sed_exposure_days +
## modSSc8: (1 | Ref)
## modSSc10: Binary_death_total ~ Sed_level_stand_mg + (1 | Gsp) + (1 | Ref)
## modSSc15: Binary_death_total ~ Sed_level_stand_mg * Sed_exposure_days +
## modSSc15: (1 | Ref_name/Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modSSc4 3 38.093 47.605 -16.0466 32.093
## modSSc8 4 13.545 26.227 -2.7726 5.545 26.548 1 2.571e-07 ***
## modSSc10 4 34.442 47.124 -13.2208 26.442 0.000 0 1
## modSSc15 6 17.545 36.568 -2.7726 5.545 20.896 2 2.900e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Data: totalSS2
## Models:
## modSSc8: Binary_death_total ~ Sed_level_stand_mg + Sed_exposure_days +
## modSSc8: (1 | Ref)
## modSSc15: Binary_death_total ~ Sed_level_stand_mg * Sed_exposure_days +
## modSSc15: (1 | Ref_name/Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modSSc8 4 13.545 26.227 -2.7726 5.5452
## modSSc15 6 17.545 36.568 -2.7726 5.5452 0 2 1
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Binary_death_total ~ Sed_level_stand_mg + Sed_exposure_days +
## (1 | Ref)
## Data: totalSS2
##
## AIC BIC logLik deviance df.resid
## 13.5 26.2 -2.8 5.5 172
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1 0 0 0 1
##
## Random effects:
## Groups Name Variance Std.Dev.
## Ref (Intercept) 5.443e-17 7.378e-09
## Number of obs: 176, groups: Ref, 8
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -862.270 967.393 -0.891 0.3728
## Sed_level_stand_mg 2.029 31.628 0.064 0.9489
## Sed_exposure_days 9.562 3.984 2.400 0.0164 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Binary_death_total ~ log10_conc + log10_dur + (1 | Ref)
## Data: totalSS2
##
## AIC BIC logLik deviance df.resid
## 13.5 26.2 -2.8 5.5 172
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1 0 0 0 1
##
## Random effects:
## Groups Name Variance Std.Dev.
## Ref (Intercept) 0 0
## Number of obs: 176, groups: Ref, 8
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2676.5 2071.5 -1.292 0.196
## log10_conc 85.2 4512.3 0.019 0.985
## log10_dur 1326.1 3518.0 0.377 0.706
## convergence code: 0
## boundary (singular) fit: see ?isSingular
##
## Call:
## glm(formula = Binary_death_total ~ Sed_level_stand_mg + Sed_exposure_days,
## family = binomial(link = "logit"), data = totalSS2)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.177 0.000 0.000 0.000 1.177
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -385.217 88625.663 -0.004 0.997
## Sed_level_stand_mg 1.064 281.353 0.004 0.997
## Sed_exposure_days 4.217 968.684 0.004 0.997
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 52.3378 on 175 degrees of freedom
## Residual deviance: 5.5452 on 173 degrees of freedom
## AIC: 11.545
##
## Number of Fisher Scoring iterations: 25
## [1] 0.9886364
##
## p 0 1
## 0 170 2
## 1 0 4
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
## Area under the curve: 0.998
## R2m R2c
## theoretical 0.9996362 0.9996362
## delta 0.9966523 0.9966523
There is no significant relationship between sediment exposure and the odds of total colony mortality (concentration: GLMM z = 0.004, p = 0.997; duration: GLMM z = 0.004, p = 0.997).
#Overlaying glmm results on figure, but prediction line is jagged...
pred_modSSc2 <- predict(modSSc2, newdata=totalSS2, type="response")
totalSS3 <- cbind(totalSS2, pred_modSSc2)
totalSS_plot <- ggplot(data = totalSS3) +
geom_point(mapping = aes(
x = Sed_level_stand_mg,
y = Binary_death_total,
color = Ref)) +
labs(x = "Sediment exposure concentration (mg/L)",
y = "Total mortality due to sediment exposure?",
color = "Study") +
scale_x_log10(limits=c(1,200),breaks=c(0.01,0.1,1,10,100,1000),
label=c("0.01","0.1","1","10","100","1000")) +
geom_line(aes(x = Sed_level_stand_mg, y = pred_modSSc2), inherit.aes=FALSE) +
theme_bw()
totalSS_plot
That plot is confusing to interpret. Let’s see if I can plot the average marginal probability, i.e., the average change in probability of the outcome across the range of the predictor of interest. This is described in some detail at the following, useful website: https://stats.idre.ucla.edu/r/dae/mixed-effects-logistic-regression/
jvalues <- with(totalSS2, seq(from = min(Sed_level_stand_mg), to = max(Sed_level_stand_mg), length.out = 100))
# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
totalSS2$Sed_level_stand_mg <- j
predict(modSSc2, newdata = totalSS2, type = "response")
})
# get the means with lower and upper quartiles
plotdat <- t(sapply(pp, function(x) {
c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))
# add in Sed_level_stand_mg values and convert to data frame
plotdat <- as.data.frame(cbind(plotdat, jvalues))
# better names and show the first few rows
colnames(plotdat) <- c("PredictedProbability", "MedianProbability",
"LowerQuantile", "UpperQuantile", "Sed_conc")
#head(plotdat)
# plot average marginal predicted probabilities
ggplot(plotdat, aes(x = Sed_conc)) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
ggplot(plotdat, aes(x = Sed_conc)) +
geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
#Overlaying average marginal predicted probabilities on figure
#by Ref
totalSS_plot2 <- ggplot() +
geom_point(data = totalSS2, mapping = aes(
x = Sed_level_stand_mg,
y = Binary_death_total,
color = Ref)) +
labs(x = "Sediment exposure concentration (mg/L)",
y = "Predicted probability of total colony\nmortality due to sediment exposure",
color = "Binary data\nby Study",
linetype = "Predicted\nProbability") +
scale_x_log10(limits=c(1,200), breaks=c(0.01,0.1,1,10,100,1000,loaelSS3),
label=c("0.01","0.1","1","10","100","1000",round(loaelSS3,digits=1))) +
geom_ribbon(data = plotdat, aes(x = Sed_conc, y = PredictedProbability,
ymin = LowerQuantile, ymax = UpperQuantile),
alpha = 0.15) +
geom_line(data = plotdat, aes(x = Sed_conc, y = PredictedProbability,
linetype = "twodash")) +
geom_line(data = plotdat, aes(x = Sed_conc, y = MedianProbability,
linetype = "solid")) +
geom_vline(xintercept=loaelSS3, linetype="dashed", color = "red") +
theme_bw() +
scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))
totalSS_plot2
jvalues <- with(totalSS2, seq(from = min(Sed_exposure_days), to = max(Sed_exposure_days), length.out = 100))
# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
totalSS2$Sed_exposure_days <- j
predict(modSSc2, newdata = totalSS2, type = "response")
})
# get the means with lower and upper quartiles
plotdat <- t(sapply(pp, function(x) {
c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))
# add in Sed_exposure_days values and convert to data frame
plotdat <- as.data.frame(cbind(plotdat, jvalues))
# better names and show the first few rows
colnames(plotdat) <- c("PredictedProbability", "MedianProbability",
"LowerQuantile", "UpperQuantile", "Sed_dur")
#head(plotdat)
# plot average marginal predicted probabilities
ggplot(plotdat, aes(x = Sed_dur)) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
ggplot(plotdat, aes(x = Sed_dur)) +
geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
#Overlaying average marginal predicted probabilities on figure
#by Ref
totalSS_plot2 <- ggplot() +
geom_point(data = totalSS2, mapping = aes(
x = Sed_exposure_days,
y = Binary_death_total,
color = Ref)) +
labs(x = expression("Sediment exposure duration (days)"),
y = "Predicted probability of total colony\nmortality due to sediment exposure",
color = "Study",
linetype = "Predicted\nProbability") +
scale_x_continuous(limits=c(0,max(totalSS2$Sed_exposure_days))) +
geom_ribbon(data = plotdat, aes(x = Sed_dur, y = PredictedProbability,
ymin = LowerQuantile, ymax = UpperQuantile),
alpha = 0.15) +
geom_line(data = plotdat, aes(x = Sed_dur, y = PredictedProbability,
linetype = "twodash")) +
geom_line(data = plotdat, aes(x = Sed_dur, y = MedianProbability,
linetype = "solid")) +
theme_bw() +
scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))
totalSS_plot2